Regresión Logística

GENERALIDADES

La regresión logística es un tipo de análisis de regresión utilizado para predecir el resultado de una variable dependiente o respuesta categórica (por lo general de dos categorías o binaria) en función de una o más variables independientes o predictoras.

La regresión logística es uno de los procedimientos estadísticos más usados en ciencias médicas, epidemiología y ciencias sociales, entre otras.

Regresión Logística

HISTORIA

La técnica de la regresión logística se originó en el año 1961 con el trabajo de Cornfield, Gordon y Smith. Pero en 1967 Walter y Duncan comenzaron a utilizarla en la forma que la conocemos actualmente, o sea para estimar la probabilidad de ocurrencia de un proceso en función de ciertas variables.

Su uso se incrementa desde principios de los 80 como consecuencia de los adelantos ocurridos en el campo de la computación

Regresión Logística

¿Por qué regresión logística y no lineal?

Si una variable cualitativa con dos niveles se codifica como 1 y 0, matemáticamente es posible ajustar un modelo de regresión lineal por mínimos cuadrados \(\beta_0+\beta_1X\) (a+bx).

El problema de esta aproximación es que, al tratarse de una recta, para valores extremos del predictor, se obtienen valores de Y menores que 0 o mayores que 1, lo que entra en contradicción con el hecho de que las probabilidades siempre están dentro del rango [0,1]

Regresión Logística

Para evitar estos problemas, la regresión logística transforma el valor devuelto por la regresión lineal \(\beta_0+\beta_1X\) empleando una función cuyo resultado está siempre comprendido entre 0 y 1.

Existen varias funciones que cumplen esta descripción, una de las más utilizadas es la función logística (también conocida como función sigmoide):

\[\huge P_(Y=k|X=x)=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}}\]

Regresión Logística

Concepto de ODDS o razón de probabilidad

En la regresión lineal simple, se modela el valor de la variable dependiente Y en función del valor de la variable independiente X.

Sin embargo, en la regresión logística, se modela la probabilidad de que la variable respuesta Y pertenezca al nivel de referencia en función del valor que adquieran los predictores, mediante el uso de LOG of ODDs.

Supóngase un evento cuya probabilidad de ser verdadero es de 0.8 (p=0,8). Dado que es una variable binaria, la probabilidad ser falso es de 1 - 0.8 = 0.2 (q=0.2).

Los odds o razón de probabilidad se definen como el ratio entre la probabilidad de evento verdadero y la probabilidad de evento falso p/q. En este caso los odds son 0.8 / 0.2 = 4, lo que equivale a decir que se esperan 4 eventos verdaderos por cada evento falso.

Regresión Logística

Coeficientes de los predictores

La estimación de los parámetros β0 y β1 se obtiene por el método de máxima verosimilitud (maximum likelihood ML), es decir, el valor de los parámetros β0 y β1 con los que se maximiza la probabilidad de obtener los datos observados.

El método de maximum likelihood está ampliamente extendido en la estadística pero su implementación y explicación no siempre es fácil de realizar y en nuestra caso particular, excede el nivel de este curso.

Afortunadamente, programas como SPSS y R, entre otros, calculan los coeficientes de forma automática.

Regresión Logística

Interpretación del modelo

Los principales elementos que hay que interpretar en un modelo de regresión logística son los coeficientes de los predictores:

β0 = es la ordenada en el origen o intercept. Se corresponde con el valor esperado del logaritmo de odds cuando todos los predictores son cero. Puede transformarse a probabilidad con la ecuación \(𝑒^{𝛽_0}/(1+𝑒^{𝛽_0} )\). Tras la transformación, su valor se corresponde con la probabilidad esperada de pertenecer a la clase 1 cuando todos los predictores son cero.

βi = los coeficientes de regresión parcial de cada predictor indican el cambio promedio del logaritmo de odds al incrementar en una unidad la variable predictora Xi, manteniéndose constantes el resto de variables. Esto equivale a decir que, por cada unidad que se incrementa Xi, se multiplican los odds por \(𝑒^{𝛽_i}\)

Regresión Logística

Condiciones para su Implementación

Independencia: las observaciones tienen que ser independientes unas de otras.

Relación lineal entre el logaritmo natural de odds y la variable continua: patrones en forma de “U” son una clara violación de esta condición.

La regresión logística no precisa de una distribución normal de la variable continua independiente.

Número de observaciones: no existe una norma establecida al respecto, pero se recomienda entre 50 a 100 observaciones. Por su parte, Freeman en 1987 ha sugerido que el número de sujetos debe ser superior a 10 * (k + 1), donde k es el número de variables independientes, o lo que es lo mismo, en términos generales, el tamaño de muestra ha de ser unas diez veces el número de variables independientes a estimar más uno.

Regresión Logística

Predicciones

Una vez estimados los coeficientes del modelo logístico, es posible conocer la probabilidad de que la variable dependiente pertenezca al nivel de referencia, dado un determinado valor del predictor. Para ello se emplea la ecuación del modelo:

\[\Large P_{(Y=1|X)}=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}}\] Convertir probabilidad en clasificación

Una de las principales aplicaciones de un modelo de regresión logística es clasificar la variable cualitativa en función de valor que tome el predictor.

Para conseguir esta clasificación, es necesario establecer un threshold (umbral) de probabilidad a partir de la cual se considera que la variable pertenece a uno de los niveles. Por ejemplo, se puede asignar una observación al grupo 1 si \(\hat p_{(Y=1|X)}>0.5\) y al grupo 0 si es lo contrario.

Regresión Logística

EJEMPLO

El fichero craneos.xlsx contiene los datos recolectados por el Coronel L.A. Waddel en el sureste y este del Tíbet y reportando por G.M. Morant en 1923.. Estos datos consisten en dos tipos de cráneos, Tipo 1 encontrados en tumbas de Sikkim y áreas vecinas del Tíbet, y los Tipo 2, recogidos en la campo de batalla de Lhasa, que se cree que fueron soldados nativos de la provincia oriental de Khans.

Se cree posible que los dos conjuntos pertenezcan a dos etnias diferentes. Cada uno de los 32 craneos tienen 5 medidas Longitud (logitud mayor del cráneo), Anchura (anchura mayor horizontal del cráneo), Altura (altura del cráneo), Altura_cara (altura de la cara superior), Anchura_cara (ancho del cara, entre los puntos más externos de los pómulos)

El link de los datos originales los puede encontrar aqui: https://math.univ-angers.fr/~labatte/enseignement%20UFR/master%20MIM/skull.pdf

Regresión Logística - Ejemplo

Cargamos los datos

library(readxl)   #Cargamos la librería
datos <- read_excel("data/craneos.xlsx")
str(datos)    #Revisamos el tipo de datos
tibble [32 x 6] (S3: tbl_df/tbl/data.frame)
 $ Longitud    : num [1:32] 19050 17250 16700 16950 17500 ...
 $ Anchura     : num [1:32] 15250 13200 13000 15050 13850 ...
 $ Altura      : num [1:32] 14500 12550 12550 13350 12600 ...
 $ Altura_Cara : num [1:32] 7350 6300 6950 6450 7750 7150 7050 7350 7000 6200 ...
 $ Anchura_Cara: num [1:32] 13650 12100 11950 12800 13550 ...
 $ Tipo_raza   : num [1:32] 1 1 1 1 1 1 1 1 1 1 ...

En seguida, tenemos que transformar la variable dependiente a 0 y 1 (asi lo requieren las regresiones logisticas binomiales)

datos$Tipo_raza <- ifelse(datos$Tipo_raza==1,0,1)
table(datos$Tipo_raza)   #verificamos el cambio
 0  1 
17 15 

Regresión Logística - Ejemplo

Para la realización de un modelo logístico, usaremos la función glm() que es una función base de R.

logist <- glm( Tipo_raza ~ Altura_Cara+Anchura_Cara, data = datos, family = binomial)
summary(logist) 

Regresión Logística - Ejemplo

Obtener los Odds ratio

Para conocer el impacto que tiene una variable independiente sobre una de la categorías de la variable dependiente en comparación con el otro, usaremos los odds ratio.

Para calcular los odds ratio en r se usa el siguiente comando

exp(coef(x))

En donde x es el objeto que contiene el modelo de la regresión logística, y exp(), calcula el valor exponencial de un número. Esto es equivlente a \(e^x\), donde e es la constante de Euler = 2.71828…

Por ejemplo:

exp(coef(logist))
 (Intercept)  Altura_Cara Anchura_Cara 
3.030043e-16 1.003117e+00 1.000951e+00 

Estos Odd, representan una medida relativa de cambio en el riesgo de un perfil respecto del otro.Para interpretarlos, consideraremos los siguiente:

  • Si los individuos comparados son idénticos en los valores que toman las variables independientes, la razón de ventajas toma el valor \(e^0=1\). Esto indica que no existen diferencias entre ellos (no hay relación).

Regresión Logística - Ejemplo

  • Valores por encima de 1, significa que la asociación es positiva entre la variable independiente y la dependiente, es decir, indica que por cada unidad que se incrementa la variable independiente, el odds de la variable dependiente (codigo 1) tambien aumenta.
  • Valores menores a 1, significa que hay una asociación negativa o inversa entre la variable independiente y la dependiente, es decir indican que que por cada unidad que se incrementa la variable independiente, el odds de la variable dependiente (codigo 1) decrece. En salud y medicina, este fenómeno se considera protector de la enfermedad o riesgo.

En este caso particular, podemos interpretar que la posibilidad que un cráneo pertenezca a la etnia Tipo 2 (codigo 1 en nuestra base), aumenta 1.0031 veces cuando se aumenta la altura de la cara en una unidad (1 mm).

Tambien podemos interpretarlo como probabilidad: (1.0031 - 1)*100 = 0.31%

paste(round((exp(coef(logist))[2]-1)*100,2), "%")
[1] "0.31 %"

Es decir, la posibilidad de pertenecer a la etnia Tipo 2 aumenta en 0.31% cuando la altura de la cara aumenta en una unidad (1 mm).

Ahora bien, si aumentamos en 10 mm (1 cms) la altura de la cara, la posibilidad de pertenecer a la etnia Tipo 2 se incrementa en: (exp(0.003112088 * 10) -1 )*100 = 3.16%

(exp(logist$coefficients[2]*10) -1)*100

Regresión Logística - Ejemplo

Obtener los Odds ratio

Para la Anchura de la cara, hacemos el mismo procedimiento, la posibilidad de que un cráneo pertenezca a la etnia Tipo 2, aumenta 1.0009 veces cuando se aumenta la anchura de la cara en una unidad (1 mm).

Tambien podemos interpretarlo como probabilidad: (1.00095 - 1)*100 = 0.095%

print(paste((exp(coef(logist))[3]-1)*100, "%"))
[1] "0.0951171037421927 %"

En este caso, vemos que la contrubución del ancho de cara es muy poca y por algo, no es una variable significativa para el modelo.

En conclusión, Altura_Cara aporta significativamente en la estimación del Tipo o afinidad racial, pero Anchura_Cara no.

Cuando los odds ratio es:

  • = 2 significa que la chance aumenta en un 100% (el doble)
  • = 3 significa que la chance aumenta en un 200% (el triple)
  • Cuando lo odds ratio son muy pequeños, menores de 1 –> hay que dividir (1/odds) p.e. si uno tiene un Odd de 0.01, entonces 1/0.01= aumenta la posibilidad en 100 veces de ser Tipo 1 (generalizando, de ser 0 )

Clasificación predicha

Como ya conocemos los coeficientes, podemos construir la ecuación regresiva y calcular posteriormente la estimación de la clasificación.

Sabemos que la ecuación de regresión es: \(\hat Y= \beta_0+\beta_1X_1+\beta_2X_2\), por tanto sustituyendo en la formula, quedaría asi: \[ \hat Y = -35.37278+(0.00311*Altura Cara)+(0.00095*Anchura Cara)\]

Comparación de clasificación estimada y observaciones

La clasificación estimada la podemos realizar con la siguiente función:

predictor <- predict(logist, newdata = datos, type = "response")
#Esta función solo nos produce la probabildad de pertenencia
head(as.data.frame(predictor)) 
    predictor
1 0.529559315
2 0.009728385
3 0.060501271
4 0.029580447
5 0.780423933
6 0.263678223

Ahora tenemos que clasificar según el umbral que nosotros, decidamos, en este caso 0.5

clasificador <- ifelse(predictor > 0.5, 1, 0)
head(as.data.frame(clasificador)) 
  clasificador
1            1
2            0
3            0
4            0
5            1
6            0

Comparación de clasificación estimada y observaciones

Ahora podemos tabular la clasificación estimada con la clasificación observada y ver las coincidencias

matriz_confusion <- table(observaciones = datos$Tipo_raza, Prediccion = clasificador)
matriz_confusion
             Prediccion
observaciones  0  1
            0 13  4
            1  4 11

Luego calculamos la Exactitud (Accuracy) como: (13+11/32)*100 = 75%. En R podemos usar este comando:

mean(clasificador==datos$Tipo_raza)*100   
[1] 75

En resumen, el modelo es capaz de clasificar correctamente el 75% de las observaciones

Finalmente realizamos la evaluación del modelo con el Test Pseudo R2 de McFadden. Es una medida de bondad de ajuste de un modelo. Se define como la proporción de variación de la variable respuesta o dependiente que es explicado por el modelo de regresión ajustado.

library(DescTools)
PseudoR2(logist)
 McFadden 
0.3703261 

GRAFICAR LA REGRESION LOGISTICA

Hay varias formas de graficas, pero la más usual es usar ggplot. En nuestro caso, como solo Altura_Cara, fué significativa, procedemos a trabajar son esa variable, ya que con dos o mas variables, es más complicado realizar el gráfico.

library(ggplot2)
ggplot(datos, aes(x=Altura_Cara, y=Tipo_raza)) + geom_point(aes(color = as.factor(Tipo_raza)), 
                                                            shape = 16, size=5)+
  geom_smooth(method = "glm", method.args = list(family = "binomial"), color = "gray20", se = FALSE) +
  theme_bw()+ theme(legend.position = "none")

ANALISIS DISCRIMINANTE

ANALISIS DISCRIMINANTE

El Análisis Discriminante Lineal o Linear Discrimiant Analysis (LDA) es un método de clasificación supervisado de variables cualitativas en el que dos o más grupos son conocidos a priori y nuevas observaciones se clasifican en uno de ellos en función de sus características. Haciendo uso del teorema de Bayes, LDA estima la probabilidad de que una observación, dado un determinado valor de los predictores, pertenezca a cada una de las clases de la variable cualitativa, P(Y=k|X=x). Finalmente se asigna la observación a la clase k para la que la probabilidad predicha es mayor.

Es una alternativa a la regresión logística cuando la variable cualitativa tiene más de dos niveles.

Cuando se trata de un problema de clasificación con solo dos niveles, ambos métodos suelen llegar a resultados similares.

Ha sido muy útil en biología para poder conocer un set de variables morfológicas que sean capaces de discriminar por ejemplo especies, dimorfismo sexual, etc.

En antropología forense han sido muy utilizadas para poder conocer el sexo o la ancestría de un individuo a partir de su morfología esqueletal.

Esencialmente, los análisis discriminantes buscan la mejor combinación de variables independientes para maximizar la distancia entre grupos.

ANALISIS DISCRIMINANTE

El proceso de un análisis discriminante puede resumirse en 6 pasos:

  • Disponer de un conjunto de datos en el que se conoce a que grupo pertenece cada observación.
  • Calcular las probabilidades previas (prior probabilities): la proporción esperada de observaciones que pertenecen a cada grupo.
  • Determinar si la varianza o matriz de covarianzas es homogénea en todos los grupos.
  • Estimar los parámetros necesarios para las funciones de probabilidad condicional, verificando que se cumplen las condiciones para hacerlo.
  • Calcular el resultado de la función discriminante. El resultado de esta determina a qué grupo se asigna cada observación.
  • Utilizar validación cruzada (cross-validation) para estimar las probabilidades de clasificaciones erróneas.

ANALISIS DISCRIMINANTE

El análisis discriminante lineal se utiliza con frecuencia como técnica de reducción de dimensionalidad para el reconocimiento o clasificación de patrones.

Los objetivos del Análisis Discriminante pueden sintetizarse en dos:

  • Analizar si existen diferencias entre los grupos en cuanto a su comportamiento con respecto a las variables consideradas y averiguar en qué sentido se dan dichas diferencias.
  • Elaborar procedimientos de clasificación sistemática de individuos de origen desconocido, en uno de los grupos analizados.

La base del análisis discriminante consiste en llegar a establecer una combinación lineal de las variables independientes que nos permita clasificar a los sujetos (poblaciones) en los diferentes grupos establecidos a priori

Esta combinación lineal es lo que se conoce como función discriminante

Hay tantas funciones discriminantes como grupos menos 1

ANALISIS DISCRIMINANTE

El origen del análisis discriminante lineal se remonta al matemático R. Fisher (1936). Su propuesta consistía en hallar la dirección e que maximizara la diferencia entre grupos respecto a la variabilidad dentro de los grupos y trabajar con estas cantidades unidimensionales.

Fisher planteó una aproximación en la que el espacio p-dimensional (donde p es el número de predictores originales) se reduce a un subespacio de menos dimensiones formado por las combinaciones lineales de las variables originales que mejor explican la separación de las clases. Una vez encontradas dichas combinaciones, se realiza la clasificación en este subespacio. El subespacio óptimo, según Fisher, es aquel que maximiza la distancia entre grupos en términos de varianza.

ANALISIS DISCRIMINANTE

Los análisis discriminantes buscan principalmente:

  1. Encontrar relaciones lineales entre las variables continuas que mejor discriminen en los grupos dados a las observaciones.

  2. Construir una regla de decisión que asigne una observción nueva con un cierto grado de riesgo, cuya clasificación previa se desconoce, a uno de los grupos prefijados (como en el caso de la antropología forense)

ANALISIS DISCRIMINANTE

Limitaciones

  • Están diseñados para variables independientes que deben distribuir normalmente, lo que lo convierte en una estrategia más limitada que las regresiones logísticas

  • No solo exigen normalidad de cada variable independiente sino también normalidad multivariante de todas las variables independientes involucradas. Esto convierte a los análisis de funciones discriminantes en una técnica muy limitada ya que es muy difícil encontrar normalidad multivariante

  • Asume igual matriz de covarianza de las variables independientes entre cada grupo

  • Requiere de un gran tamaño muestral para obtener coeficientes insesgados (no ocurre así con la regresión logística)

ANALISIS DISCRIMINANTE

Validación del LDA

Contraste de normalidad Shapiro-Wilk o kolmogorov-smirnov para cada variable

Contraste de normalidad multivariante de Henze-Zirkler, Doornik-Hansen o Royston con el paquete El paquete MVN

La homogeneidad de varianzas se puede realizar con el El test Box M (1949) del paquete biotools. El test Box M es muy sensible a violaciones de la normalidad multivariante, por lo que se recomienda emplear un límite de significancia de 0.001

Matriz de Confusión para comparar la clasificación estimada con la clasificación observada. El porcentaje de individuos correctamente clasificados nos va a dar información sobre la confianza de la ecuación predictiva.

Generalmente en Antropología Forense se considerará una ecuación confiable aquella que tenga al menos un 85% de casos correctamente clasificados.

ANALISIS DISCRIMINANTE

EJEMPLO

Usaremos el mismo fichero craneos.xlsx con los datos del Tibet

library(readxl)
craneos <- read_excel("data/craneos.xlsx")   #Ahora le llamamos craneos al objeto
#Revisamos la data
str(craneos)
tibble [32 x 6] (S3: tbl_df/tbl/data.frame)
 $ Longitud    : num [1:32] 19050 17250 16700 16950 17500 ...
 $ Anchura     : num [1:32] 15250 13200 13000 15050 13850 ...
 $ Altura      : num [1:32] 14500 12550 12550 13350 12600 ...
 $ Altura_Cara : num [1:32] 7350 6300 6950 6450 7750 7150 7050 7350 7000 6200 ...
 $ Anchura_Cara: num [1:32] 13650 12100 11950 12800 13550 ...
 $ Tipo_raza   : num [1:32] 1 1 1 1 1 1 1 1 1 1 ...
summary(craneos[ , c(1:5)])
    Longitud        Anchura          Altura       Altura_Cara    Anchura_Cara  
 Min.   :16250   Min.   :12650   Min.   :12150   Min.   :6200   Min.   :11850  
 1st Qu.:17325   1st Qu.:13488   1st Qu.:12962   1st Qu.:6988   1st Qu.:13138  
 Median :17950   Median :13900   Median :13300   Median :7350   Median :13450  
 Mean   :17994   Mean   :13906   Mean   :13330   Mean   :7294   Mean   :13370  
 3rd Qu.:18463   3rd Qu.:14250   3rd Qu.:13662   3rd Qu.:7650   3rd Qu.:13662  
 Max.   :20000   Max.   :15300   Max.   :14500   Max.   :8250   Max.   :14650  

ANALISIS DISCRIMINANTE EJEMPLO

Convertimos factor la variable Tipo_raza para asegurarnos que sea leida como categórica

craneos$Tipo_raza <- factor(craneos$Tipo_raza)

Realizamos los graficos QQ-plot para probar normalidad

par(mfrow=c(2,2))
for(i in 1:4){
  qqnorm(unlist(craneos[,i]), pch = 1)
  qqline(unlist(craneos[,i]), col = "steelblue", lwd = 2)}

par(mfrow=c(1,1))

ANALISIS DISCRIMINANTE EJEMPLO

Visualizamos el último gráfico

qqnorm(unlist(craneos[,5]), pch = 1)
qqline(unlist(craneos[,5]), col = "steelblue", lwd = 2)

ANALISIS DISCRIMINANTE EJEMPLO

Ahora realizamos los Test de Normalidad para cada variable

by(craneos$Longitud, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.97442, p-value = 0.8901

------------------------------------------------------------ 
craneos$Tipo_raza: 2

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.94141, p-value = 0.4004
by(craneos$Anchura, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.96081, p-value = 0.6468

------------------------------------------------------------ 
craneos$Tipo_raza: 2

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.94619, p-value = 0.4666

ANALISIS DISCRIMINANTE EJEMPLO

by(craneos$Altura, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.9586, p-value = 0.6053

------------------------------------------------------------ 
craneos$Tipo_raza: 2

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.9654, p-value = 0.785
by(craneos$Altura_Cara, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.96785, p-value = 0.7798

------------------------------------------------------------ 
craneos$Tipo_raza: 2

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.97646, p-value = 0.9396

ANALISIS DISCRIMINANTE EJEMPLO

by(craneos$Anchura_Cara, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.94674, p-value = 0.4071

------------------------------------------------------------ 
craneos$Tipo_raza: 2

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.93485, p-value = 0.322

Ahora realizamos el Test de Normalidad Multiple con el paquete MVN

library(MVN)
mvn(craneos[ , c(1:5)],mvnTest = "hz")  #Test de Henze-Zirkler
$multivariateNormality
           Test        HZ   p value MVN
1 Henze-Zirkler 0.8261507 0.3755845 YES

$univariateNormality
              Test     Variable Statistic   p value Normality
1 Anderson-Darling   Longitud      0.4211    0.3048    YES   
2 Anderson-Darling   Anchura       0.3057    0.5470    YES   
3 Anderson-Darling    Altura       0.2250    0.8049    YES   
4 Anderson-Darling Altura_Cara     0.1872    0.8962    YES   
5 Anderson-Darling Anchura_Cara    0.7763    0.0392    NO    

$Descriptives
              n     Mean  Std.Dev Median   Min   Max    25th    75th       Skew
Longitud     32 17993.75 936.5129  17950 16250 20000 17325.0 18462.5  0.4267349
Anchura      32 13906.25 684.1229  13900 12650 15300 13487.5 14250.0  0.2908175
Altura       32 13329.69 608.2576  13300 12150 14500 12962.5 13662.5  0.0930111
Altura_Cara  32  7293.75 539.0778   7350  6200  8250  6987.5  7650.0 -0.2261243
Anchura_Cara 32 13370.31 744.4265  13450 11850 14650 13137.5 13662.5 -0.4722848
               Kurtosis
Longitud     -0.6532942
Anchura      -0.6738678
Altura       -0.8202897
Altura_Cara  -0.8175725
Anchura_Cara -0.3977022

ANALISIS DISCRIMINANTE EJEMPLO

Aplicamos la función boxM() del paquete biotools que realiza el contraste de homogeneidad de matrices de covarianza.

#devtools::install_github("arsilva87/biotools")  #Para instalar este paquete usamos este código
library(biotools)

boxM(data = craneos[, 1:5], grouping = craneos$Tipo_raza)
    Box's M-test for Homogeneity of Covariance Matrices

data:  craneos[, 1:5]
Chi-Sq (approx.) = 18.371, df = 15, p-value = 0.2437

Observamos que no tenemos elementos para rechazar la homogeneidad de covarianzas

ANALISIS DISCRIMINANTE EJEMPLO

Ahora que todos los test nos dan bien, procedemos a realizar el Analisis Discriminante Lineal con el paquete MASS

#install.packages("MASS")     #Recuerden instalarlo primero
library(MASS)
modelo_lda <- lda(Tipo_raza ~ Longitud+Anchura+Altura+Altura_Cara+Anchura_Cara, data = craneos)

ANALISIS DISCRIMINANTE EJEMPLO

Podemos graficar el modelo

plot(modelo_lda)

ANALISIS DISCRIMINANTE EJEMPLO

Con los coefiencientes de la función discriminante, construimos la ecuación predictiva:

\(\hat Y=(0.0004772 * Longitud)+(-0.0008324*Anchura)+(-0.0000279*Altura)+\) \((0.0009469*\verb|Altura_Cara|)+(0.000948*\verb|Anchura_Cara|)\)

Solo nos falta la constante, que el paquete MASS no lo muestra automáticamente, aunque internamente si lo usa para los calculos. Nosotros podemos obtenerlo con el siguiente código:

groupmean<-(modelo_lda$prior %*% modelo_lda$means)  #operador % * % es para realizar producto de dos matrices.
constant<--(groupmean %*% modelo_lda$scaling)
constant
           LD1
[1,] -16.22159

Para ayudarnos a clasificar los valores o scores de la ecuación predictiva, usamos como referencia los centroides de las funciones. Para eso, sacamos los promedios de los scores predictivos

lda_prediction <- predict(modelo_lda)   #Aqui obtenemos los scores y clases pronosticadas

tapply(lda_prediction$x, craneos$Tipo_raza, mean)  #Así obtenemos las funciones de los centroides de los grupos
         1          2 
-0.8771315  0.9940824 

ANALISIS DISCRIMINANTE EJEMPLO

Los cálculos de los scores y clase pronosticada se realizan de la siguiente forma:

Todo esto se guarda en el objeto “lda_prediction”, y tiene 3 elementos importantes:

  • lda_prediction$x –> posee los scores pronosticados
  • lda_prediction$posterior –> Posee las probabilidades de pertenencia para cada grupo
  • lda_prediction$class –> La clase pronosticada. Usaremos este para evaluar la exactitud del modelo

ANALISIS DISCRIMINANTE EJEMPLO

Para la exactitud del modelo, usaremos una tabla de contingencia o matriz de confusión para tabular la clasificación estimada con la clasificación observada y ver las coincidencias

table(predict=lda_prediction$class, Actual=craneos$Tipo_raza)
       Actual
predict  1  2
      1 14  3
      2  3 12

Luego calculamos la Exactitud (Accuracy) como: (14+12/32)*100 = 81.25%.

mean(lda_prediction$class==craneos$Tipo_raza)*100
[1] 81.25

En resumen, el modelo clasificó correctamente el 81.25% de las observaciones

ANALISIS DISCRIMINANTE EJEMPLO

Graficos de Partición por cada par de variables

#install.packages("klaR")  #Recuerden instalarlo primero
library(klaR)
partimat(craneos[, 1:5],grouping=as.factor(craneos$Tipo_raza), method = "lda")

ANALISIS DISCRIMINANTE EJEMPLO

Finalmente, graficamos, como en este caso tenemos una sola función discriminante, procedemos así:

prediction <- data.frame(lda_prediction)   #Guardamos nuestro objeto de predicción como un data frame

ggplot(prediction, aes(x=as.factor(class) , y=LD1, fill = as.factor(class)))+
  geom_dotplot(binaxis = "y", stackdir = "center", binwidth = 0.25)+
  xlab("Pronostico") + theme_light() + theme(legend.position = "none")

ANALISIS DISCRIMINANTE EJEMPLO

EJEMPLO CON MÁS DE DOS GRUPOS

Para este ejemplo, tomaremos los datos de las frecuencias sanguíneas de ABO, disponibles en la página de Wikipedia

Posteriormente, se limpió y se agregó los continentes de cada país.

Blood_type_new <- read_excel("data/Blood_type_final.xlsx")
head(Blood_type_new)
# A tibble: 6 x 11
  continent country   population  `O+`  `A+`  `B+`  `AB+`  `O-`   `A-`  `B-`
  <chr>     <chr>          <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl>
1 Europa    Albania      3074579 0.341 0.312 0.145 0.052  0.06  0.055  0.026
2 Africa    Algeria     43576691 0.4   0.3   0.15  0.0425 0.066 0.023  0.011
3 America   Argentina   45479118 0.489 0.315 0.08  0.0245 0.049 0.0316 0.008
4 Asia      Armenia      3021324 0.29  0.463 0.12  0.056  0.02  0.037  0.01 
5 Oceania   Australia   25466459 0.38  0.32  0.12  0.04   0.07  0.06   0.02 
6 Europa    Austria      8859449 0.3   0.37  0.12  0.05   0.06  0.07   0.02 
# i 1 more variable: `AB-` <dbl>

Nos aseguramos que la variable continente sea leida como categórica

Blood_type_new$continent <- factor(Blood_type_new$continent)

Y eliminamos a Oceanía por tener poca frecuencia

Blood_type_new <- Blood_type_new[Blood_type_new$continent != "Oceania", ] 

ANALISIS DISCRIMINANTE EJEMPLO

Revisamos las variables para asegurarnos que no existan valores perdidos o extraños

summary(Blood_type_new[ , c(4:11)])
       O+               A+               B+              AB+         
 Min.   :0.2680   Min.   :0.1400   Min.   :0.0472   Min.   :0.00500  
 1st Qu.:0.3247   1st Qu.:0.2479   1st Qu.:0.1000   1st Qu.:0.02900  
 Median :0.3817   Median :0.2990   Median :0.1510   Median :0.04295  
 Mean   :0.4044   Mean   :0.2956   Mean   :0.1637   Mean   :0.04756  
 3rd Qu.:0.4684   3rd Qu.:0.3507   3rd Qu.:0.2114   3rd Qu.:0.06313  
 Max.   :0.7500   Max.   :0.4630   Max.   :0.3680   Max.   :0.14400  
       O-                A-                B-               AB-          
 Min.   :0.00060   Min.   :0.00040   Min.   :0.00010   Min.   :0.000100  
 1st Qu.:0.01798   1st Qu.:0.01000   1st Qu.:0.00585   1st Qu.:0.001225  
 Median :0.04000   Median :0.02715   Median :0.01275   Median :0.004150  
 Mean   :0.03907   Mean   :0.03284   Mean   :0.01340   Mean   :0.004937  
 3rd Qu.:0.06000   3rd Qu.:0.06000   3rd Qu.:0.02000   3rd Qu.:0.009025  
 Max.   :0.13000   Max.   :0.08000   Max.   :0.03130   Max.   :0.012000  

Para trabajar menos variables, podemos sumar los grupos O, los grupos A, los grupos B y finalmente los AB, es decir, no tomar en cuenta los Rh.

Blood_type_new$O <- Blood_type_new$`O+`+Blood_type_new$`O−`
Blood_type_new$A <- Blood_type_new$`A+`+Blood_type_new$`A−`
Blood_type_new$B <- Blood_type_new$`B+`+Blood_type_new$`B−`
Blood_type_new$AB <- Blood_type_new$`AB+`+Blood_type_new$`AB−`

Para efectos practicos vamos obviar las pruebas de normalidad. Pero quedan de tarea.

ANALISIS DISCRIMINANTE EJEMPLO

Realizamos el análisis discriminante

modelo_lda <- lda(continent ~ O+A+B+AB, data = Blood_type_new)
Warning in lda.default(x, grouping, ...): group Oceania is empty
modelo_lda
Call:
lda(continent ~ O + A + B + AB, data = Blood_type_new)

Prior probabilities of groups:
   Africa   America      Asia    Europa 
0.1750000 0.1583333 0.3500000 0.3166667 

Group means:
                O         A         B         AB
Africa  0.5032857 0.2672524 0.1932762 0.04118571
America 0.5735368 0.2982684 0.1036421 0.02454737
Asia    0.3974405 0.2963714 0.2395762 0.06745714
Europa  0.3963421 0.4127105 0.1356842 0.05617105

Coefficients of linear discriminants:
        LD1      LD2      LD3
O  16.99326 30.00720 65.17075
A  14.39362 12.43382 74.58440
B  33.72184 19.74911 78.34340
AB 22.18339 21.79194 24.73584

Proportion of trace:
   LD1    LD2    LD3 
0.5183 0.4659 0.0158 

ANALISIS DISCRIMINANTE EJEMPLO

Como tenemos más de dos funciones discriminantes, podemos graficar las dos primeras:

plot(modelo_lda, dimen=2, col = as.integer(Blood_type_new$continent))

ANALISIS DISCRIMINANTE EJEMPLO

Verificamos la exactitud del modelo

lda_prediction <- predict(modelo_lda)

table(predict=lda_prediction$class, Actual=Blood_type_new$continent)
         Actual
predict   Africa America Asia Europa Oceania
  Africa       8       2    5      0       0
  America      3      13    0      2       0
  Asia         8       0   28      3       0
  Europa       2       4    9     33       0
  Oceania      0       0    0      0       0

Entonces la Exactitud (Accuracy) es: (8+13+28+33/120)*100 = 68.3%%.

mean(lda_prediction$class==Blood_type_new$continent)*100
[1] 68.33333

El modelo clasifica correctamente el 68.3% de las observaciones